Code:
OUTPUT: for reading in entry #3 through #68 when n = 1 (2nd time) ---
$eqparms_42000 19.801934 19.987543 19.998929 16.280022
15.882664 11.936946 11.591798 -0.069747 -0.115387 -0.036517
-0.056792 -0.040874 -0.020368 0.018431 0.014045 0.006112
0.016178 -0.160537 0.010308 0.033332 -0.012482 0.013886
-0.075848 -0.016654 -0.076822 0.028593 -0.022062 -0.109674
0.011751 0.004857 -0.024652 0.035485 -0.029861 -0.019982
-0.091736 0.036087 -0.067619 -0.018104 -0.138921 -0.109324
-0.103670 -0.140229 0.011544 -0.104651 0.071105 0.127234
0.175506 0.102106 -0.005422 13.531374 -3.500000 -3.500000
-3.500000 -3.500000 8.158419 0.199930 0.156771 0.199964
0.199996 0.199999 0.130350 14.264143 6.867160 8.385001
9.316451 14.153305 13.028459 [ it's OK here.]
/**** NOT OK at all for A[0] (long number) and A[1] (underlined) ****/
1 0 10294483105380924209297137052132835902603774975587521311686
85968221259637419742397639086525069280384040936805581819179
12296253896137979821024575318480801440811370149050155305270
13827695845672837000807447974224711286272722947708347189945
6206847082496.000000 19.987543
1 1 2.023529 19.998929
=======
/**** NOT OK at all *****/
these values of A[i] are fine.
++++++++
1 2 16.280637 16.280022
1 3 15.881519 15.882664
1 4 11.940799 11.936946
1 5 11.598596 11.591798
1 6 -0.069730 -0.069747
1 7 -0.115599 -0.115387
1 8 -0.036446 -0.036517
INPUT:
$eqparms_32000 19.809179201986 19.97629013 19.99800900 16.28063724 15.88151892 11.94079894 11.59859621 -0.06972984 -0.11
559869 -0.03644595 -0.05680799 -0.04092679 -0.02044219 0.01844713 0.01405800 0.00612413 0.01619741 -0.16061326
0.01028455 0.03333264 -0.01248030 0.01388489 -0.07584439 -0.01666025 -0.07679878 0.02859141 -0.02205941 -0.10967
226 0.01174919 0.00484831 -0.02461987 0.03549759 -0.02984242 -0.02000749 -0.09174563 0.03610507 -0.06772237 -0
.01811469 -0.13900662 -0.10938244 -0.10355774 -0.13992939 0.01139550 -0.10472170 0.07108921 0.12723564 0.17551409
0.10207980 -0.00539671 13.53033691 -3.50000000 -3.50000000 -3.50000000 -3.50000000 8.15659611 0.19995521 0.15
501831 0.19998754 0.19997583 0.19996031 0.12685206 14.27656879 6.86580210 8.38285901 9.31949826 14.15059040 1
3.03651790
$eqparms_42000 19.801934046965 19.98754311 19.99892857 16.28002229 15.88266435 11.93694552 11.59179797 -0.06974732 -0.11
538692 -0.03651734 -0.05679196 -0.04087373 -0.02036764 0.01843055 0.01404479 0.00611223 0.01617807 -0.16053726
0.01030812 0.03333178 -0.01248197 0.01388638 -0.07584773 -0.01665407 -0.07682184 0.02859282 -0.02206247 -0.10967
386 0.01175095 0.00485716 -0.02465235 0.03548512 -0.02986135 -0.01998163 -0.09173565 0.03608734 -0.06761861 -0
.01810435 -0.13892145 -0.10932366 -0.10367009 -0.14022905 0.01154401 -0.10465116 0.07110506 0.12723429 0.17550641
0.10210584 -0.00542154 13.53137379 -3.50000000 -3.50000000 -3.50000000 -3.50000000 8.15841883 0.19992979 0.15
677102 0.19996408 0.19999581 0.19999900 0.13035006 14.26414319 6.86715977 8.38500113 9.31645129 14.15330546 1
3.02845903
/* dimCov.c
**
** Read data from .ONE file, do their matrix sum,
** compute covariance and correlation.
*/
#include <math.h>
#include <stdio.h>
#include <limits.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h> /* for getopt */
#include "error.h"
#define SIZE 66
/* ---- Function Prototypes ---- */
/* Read the title for the eqparms. */
char *read_title(FILE *fp);
/* Read the rms for the eqparms. */
double read_index(FILE *fp);
/* Read entries from a line in given datafile. */
double *read_entries(FILE *fp,int fix);
double *AddMatrix(double *A, double *B);
/* compute covariance of vectors i and j */
double covar(int i,int j);
/* Global variables */
double **mat, cov[15][15], correl[15][15];
double *A; /* holds the sum and used for avg */
/*** Constants *************************************************/
#define OPTS ":cChrR" /* command line option string */
/******* Help and usage ***************************************/
static const char usage[] =
"Usage: dcov [-c][-C][-h][-r][-R] <datafile>\n";
static const char help[] =
"Usage: dcov [options] <datafile>\n\n"
"Arguments:\n"
" <datafile> datafile for which we evaluate the covariance\n"
"\t\t\tmatrices of the eqparms\n"
"Options:\n"
" -c prints covariance matrix in readable form\n"
" -C prints covariance matrix as a line (default)\n"
" -h prints this help message\n"
" -r prints correlation matrix in readable form\n"
" -R prints correlation matrix as a line (default)\n\n"
"Please report bugs to <[email protected]>. Thank you!\n";
/*** Main program ***********************************************/
int main(int argc, char **argv)
{
int c; /* used to parse command line options */
char *filename; /* name of input datatfile */
char *tempfile; /* name of temporary file */
char *shell_cmd; /* command line pointer */
char *til;
FILE *fp; /* pointer to the datafile */
int Nmax;
int i,n; /* counters */
double *rms; /* pointer to rms values */
char **Name; /* pointer to pointers of strings */
/* the NEXT few variables are read as arguments to command line options */
int lcovar = 1; /* flag for printing cov matrix as a line */
int r_LIN = 1; /* flag for printing corr_matrix as a line */
/* external declarations for command line option parsing (unistd.h) */
extern char *optarg; /* command line option argument */
extern int optind; /* pointer to current element of argv */
extern int optopt; /* contain option character upon error */
/* following part parses command line for options and their arguments */
optarg = NULL;
while ( (c = getopt(argc, argv, OPTS)) != -1 )
switch (c) {
case 'c': /* prints covariance matrix in readable form */
lcovar = 0;
break;
case 'C': /* prints covariance matrix as a line */
lcovar = 1;
break;
case 'h': /* -h help option */
PrintMsg(help, 0);
break;
case 'r':
r_LIN = 0; /* prints correlation matrix in readable form */
break;
case 'R': /* prints correlation matrix as a line */
r_LIN = 1;
break;
case ':':
error("printscore: need an argument for option -%c", optopt);
break;
case '?':
default:
error("rone: unrecognized option -%c", optopt);
}
/* error check */
if ( (argc - (optind - 1 )) != 2 )
PrintMsg(usage, 1);
/*************************************************************************
** Some file processing stuff.
*************************************************************************
*************************************************************************/
filename = (char *)calloc(CHAR_MAX, sizeof(char));
tempfile = (char *)calloc(CHAR_MAX, sizeof(char));
filename = argv[optind];
sprintf(tempfile,"%s.%s", filename,"temp");
/* Use word count and know how much data to handle.*/
shell_cmd = (char *)calloc(CHAR_MAX, sizeof(char));
sprintf(shell_cmd,"wc -l %s > %s", filename, tempfile);
if ( -1 == system(shell_cmd) ){
printf("error with shell_cmd %s", shell_cmd);
exit(1);
}
/* first, determine how many matrices there are to be stored! */
fp = fopen(tempfile,"r");
/* reads the value of n and the name of the temp file*/
til = (char *)calloc(CHAR_MAX, sizeof(char));
fscanf(fp,"%d",&Nmax);
fscanf(fp,"%s",til);
fclose(fp);
printf(" '%s' has %d lines. \n\n",filename,Nmax);
/* remove temporary file once it is no longer needed */
sprintf(shell_cmd,"rm -f %s ", tempfile);
if ( -1 == system(shell_cmd) ){
printf("error with removing temp file: %s\n", shell_cmd);
exit(1);
}
/****************************************************************
** Read the matrices from datafile
****************************************************************
****************************************************************/
fp = fopen(filename,"r"); /* Get space as needed for variables */
Name = (char **)calloc(Nmax, sizeof(char *));
mat = (double **)calloc(Nmax, sizeof(double *));
rms = (double *)calloc(Nmax, sizeof(double));
A = (double *)calloc(SIZE, sizeof(double )); /* Do sum then average with A */
for (n=0;n<Nmax;n++){
Name[n] = (char *)calloc(CHAR_MAX, sizeof(char));
til = read_title(fp);
Name[n] = strcpy(Name[n],til);
printf("%s\t",Name[n]);
rms[n] = read_index(fp);
printf("%lf ",rms[n]);
mat[n] = (double *)calloc(SIZE, sizeof(double));
mat[n] = read_entries(fp,n);
for (i=0;i<SIZE;i++)
printf("% 2.6lf ",mat[n][i]);
printf("\n\n");
for (i=0;i<SIZE;++i){
printf("%2d %2d % 2.6lf % 2.6lf\n",n,i, A[i], mat[n][i]);
A[i] += mat[n][i];}
printf("\n");
}
printf("\n");
/* USE LATER for covariance.
mat = (double **)calloc(Nmax,sizeof(double *));
for (n=0;n<Nmax;n++)
mat[n] = (double *)calloc(SIZE, sizeof(double ));
*/
if ( lcovar ){
printf("\t Printing matrix as a one-liner...\n");
printf("$Matrix_SUM A:\t");
/* for (i=0; i<SIZE;i++)
printf("% 2.6lf ", A[i]);
printf("\n\n"); */
}
else
printf("Covariance matrices printed in readable form.\n");
if ( r_LIN )
printf("Printing A as a one-liner & for each matx M, M - A...\n");
else
printf("Correlation matrices printed in readable form.\n");
/****************************************************************
** Compute the matrice A whose entries are the averages
** of homologous entries of matrices read from the datafile.
****************************************************************
****************************************************************/
/*
for (i=0;i<Nmax;i++)
A[i] /= Nmax;
printf("$Matrix_Avg A:");
for (i=0; i<SIZE;i++){
A[i] /= Nmax; take average
printf(" % 2.6lf ", A[i]);
}
printf("\n\n");
*/
/****************************************************************
** For each matrix M_i, Subtract M_i - A
****************************************************************
****************************************************************/
/*
for (n=0; n < Nmax;n++){
for (i=0; i<15;i++){
mat[n][i] -= A[i]; subtract
printf("% 1.2lf ", mat[n][i]);}
printf("\n");
}
printf("\n");
*/
/****************************************************************
** Compute Covariance of the column vectors in the 3 matrices
****************************************************************
****************************************************************/
/*
printf("Covariance matrix: Upper Half Triangle Printed .\n");
for (i=0; i<15;i++) initialize
for (j=0; j<15;j++)
{cov[i][j]=0; correl[i][j]=0;}
for (i=0; i<15;i++){
for (j=0; j<15;j++){
if ( j < i)
printf("\t");
else{ compute covariance
cov[i][j] = covar(i,j);
printf("% *.*lf\t",1,3,cov[i][j]);
}
}
printf("\n");
}
printf("\n");
*/
/****************************************************************
** Compute Covariance of the column vectors in the 3 matrices
****************************************************************
****************************************************************/
/*
printf("Correlation matrix: Upper Half Triangle Printed .\n");
for (i=0; i<15;i++){
for (j=0; j<15;j++){
if ( j < i)
printf("\t");
else{
correl[i][j] = covar(i,j)/(covar(i,i)*covar(j,j));
printf("% *.*lf\t",3,1,correl[i][j]);
}
}
printf("\n");
}
*/
/* Clean Up */
fclose(fp);
free(rms);
free(Name);
free(filename); /* name of input datatfile */
free(tempfile); /* name of temporary file */
free(shell_cmd); /* command line pointer */
free(til);
free(optarg);
return 0;
}
/* Read the title for the eqparms. */
char *read_title(FILE *fp)
{
char *title;
fscanf(fp,"%s \t",title);
/* printf("title %s \t", title);*/
return title;
}
/* Read the rms for the eqparms. */
double read_index(FILE *fp)
{
double indx;
fscanf(fp,"%lf", &indx);
/* printf("index %lf\n", indx);*/
return indx;
}
/* Read entries from a line in given datafile. */
double *read_entries(FILE *fp,int fix)
{
int i;
double *B;
B = (double *)calloc(SIZE, sizeof(double));
for (i=0;i<SIZE;i++)
fscanf(fp,"%lf", &B[i]);
return B;
}
double *AddMatrix(double *A, double *B)
{ int i;
double D[SIZE];
for (i=0;i<SIZE;i++)
D[i]=A[i]+B[i];
return D;
}
/*******************************************************
** Computes the covariances of column vectors i and j *
** Compute Covariance of the column vectors in the *
** 3 matrices *
*******************************************************
*******************************************************/
double covar(int i,int j)
{ int k;
double covaria=0.0;
for (k=0;k<3;k++) /* compute covariance */
covaria += (mat[k][i]*mat[k][j])/3;
return covaria;
}
